Predicting River Resilience to Disturbance Using Population Models at a Continental Scale

Purpose Statement

We are interested in learning about streams’ responses to disturbance (i.e., storm events). Using a state-space time series approach, we are fitting a logistic population growth model to gross primary production estimates to gain information about growth parameters (i.e., maximum growth rate, critical disturbance thresholds) that describe algal population dynamics and compare results across stream sites.

Updated Hypotheses

Q1: How does maximum growth rate (rmax) vary in streams and rivers across the country as a function of flow disturbance and light availability?

H1A: The maximum growth rate will be greatest in streams with fewer flow disturbance events that scour existing biomass, since gross primary production (GPP) tends to be highest in streams with more predictable flow (Bernhardt et al., 2022).

H1B: As light availability increases, the maximum growth rate will increase, since GPP tends to be highest in well-lit streams (Bernhardt et al., 2022).

Q2: How does the critical discharge threshold at which biomass is removed (c) vary in streams and rivers across the country as a function of flow disturbance and dominant watershed land use?

H2A: The critical discharge threshold will be greatest in streams with fewer flow disturbance events that scour existing biomass, since these sites will have likely higher maximum growth rates (see H1A) and higher GPP overall.

H2B: As land use changes relative to undeveloped watersheds, the critical discharge threshold will change according to how the dominant land use in a watershed changes flow regimes. In predominantly forested or undeveloped watersheds, critical discharge thresholds will be greater, because flow will be less flashy and more predictable (Poff et al., 2006). In predominantly urban watersheds, critical discharge thresholds will be lower, because flow will be flashier and less predictable, in large part due to higher impervious land cover (Blaszczak et al., 2019).

Model Structure

Since our last meeting, we experimented with structuring the model in a hierarchical format so as to share information across sites. However, this approach ended up being too computationally intensive (i.e., would have taken >1 month to run), so we have reverted to modeling each site individually using the original model formula (Figure 1).

Figure 1. The state-space population model being used to estimate parameters of maximum growth rate (rmax) across individual stream sites. The observation model uses daily estimates of gross primary production (gt) and daily light (Lt) to estimate daily biomass (bt). The daily persistence term (Pt) is a function of daily discharge (Qt) which is used to estimate the sensitivity of the persistence transition from presence to removal of biomass (s) and the critical discharge threshold at which benthic biomass is removed (c). Both the observation and persistence models then inform the process model, which is a density-dependent growth function for daily biomass (bt) in a given stream. The estimated terms of Pt and bt are input into the process model to extract values for the maximum growth rate (rmax) and lambda(λ), which is equal to the negative maximum growth rate divided by the carrying capacity (K).

Data Acquisition

Data have been assembled using the datasets available in Appling et al., 2018, Blaszczak et al., 2021, Bernhardt et al., 2022, and Savoy and Harvey, 2021.

We have filtered the available data according to the following rules:

  • Filter for high-quality sites based on model diagnostics (e.g., K600 sigma \(\hat{R}\) < 1.05, GPP < 0 or ER > 0 for more than 15% of estimates).
  • Filter for sites that have a maximum time gap of 14 days and a minimum of 90 days per year of daily GPP estimates.
  • Any GPP estimates below zero were assigned a value between 0.05 and 0.13.

From the initial 356 river sites in the Appling dataset, our filtered dataset consists of 182 sites and 691 site-years of data (Table 1).

Table 1. Summary statistics for all 182 stream sites included in the dataset.

Table 1. Summary statistics for all 182 stream sites included in the dataset.

Results

Figures of all available covariates (discharge, temperature, light) as well as the full set of our model diagnostics are available to view on a separate R Shiny application located here. Below are a few overview figures examining relationships among model parameter values and site covariates as they pertain to the hypotheses above. Additional site covariates are presented in the supplementary figures below.

Prior to plotting, all sites whose maximum growth rate parameter (rmax) estimates were negative (i.e., not biologically feasible) or rmax had an \(\hat{R}\) value greater than 1.05 have been removed (n = 159 sites remaining).

Main Takeaway #1: Maximum growth rate (rmax) values tend to be greater in less frequently disturbed streams (i.e., with lower coefficient of variation in discharge values, Fig. 2C).

Figure 2. (A) Distribution of median maximum growth rate (r~max~) estimates as well as plots of median r~max~ values versus (B) mean daily gross primary production (GPP), (C) the coefficient of variation in discharge (CV~Q~), (D) cumulative summer light availability, (E) road density in the catchment, and (F) percent impervious land cover in the catchment. Note, panel B presents the x-axis on a log scale.

Figure 2. (A) Distribution of median maximum growth rate (rmax) estimates as well as plots of median rmax values versus (B) mean daily gross primary production (GPP), (C) the coefficient of variation in discharge (CVQ), (D) cumulative summer light availability, (E) road density in the catchment, and (F) percent impervious land cover in the catchment. Note, panel B presents the x-axis on a log scale.

Typically, the model performs well when predicting declines in GPP (due to its design), but it underpredicts the rate at which these systems recover following disturbance (Figure 3).

Figure 3. A comparison of gross primary productivity (GPP) estimates used during the initial model fit (dark green points) versus re-predicted GPP values using our model parameter estimates (light green line with 2.5% and 97.5% confidence intervals shaded). Note, this represents only one site (N.E. Anacostia River in Riverdale, MD) for one year (2011).

Figure 3. A comparison of gross primary productivity (GPP) estimates used during the initial model fit (dark green points) versus re-predicted GPP values using our model parameter estimates (light green line with 2.5% and 97.5% confidence intervals shaded). Note, this represents only one site (N.E. Anacostia River in Riverdale, MD) for one year (2011).

Prior to plotting, all sites whose critical disturbance threshold (c) estimates were negative (i.e., not hydrologically feasible) or c had an \(\hat{R}\) value greater than 1.05 have been removed (n = 141 sites remaining).

Main Takeaway #2: In most streams, the discharge threshold at which biomass is disturbed (Qc) occurs far more frequently than the two-year flood (Q2yr) (Fig. 4).

Figure 4. Plots of the ratio of the median critical discharge (Q~c~) to the two-year flood (Q~2yr~) as (A) an overall distribution across all sites and versus (B) mean daily gross primary production (GPP), (C) the coefficient of variation in discharge (CV~Q~), (D) watershed area, (E) road density in the catchment, and (F) percent impervious land cover in the catchment. Note, the horizontal black line indicates when the Q~c~:Q~2yr~ ratio value is one; values below this line indicate the critical discharge threshold at a given site falls below the 2-year flood discharge value. Also, panels B and D present the x-axis on a log scale.

Figure 4. Plots of the ratio of the median critical discharge (Qc) to the two-year flood (Q2yr) as (A) an overall distribution across all sites and versus (B) mean daily gross primary production (GPP), (C) the coefficient of variation in discharge (CVQ), (D) watershed area, (E) road density in the catchment, and (F) percent impervious land cover in the catchment. Note, the horizontal black line indicates when the Qc:Q2yr ratio value is one; values below this line indicate the critical discharge threshold at a given site falls below the 2-year flood discharge value. Also, panels B and D present the x-axis on a log scale.

Next Steps: Following this meeting, we intend to pursue two main activities:

  1. Finalize model diagnostics to determine the list of sites included in the results/discussion.

  2. Finish a first draft of the manuscript by the end of the calendar year. Please see the working Google document here.

Questions for You:

  1. Do you have any suggestions regarding model diagnostic best practices?

  2. How should we caveat sites that are on the same river? If they are assigned different USGS site IDs, is that sufficient to say they are encountering different enough environmental conditions?

  3. Do you feel there are particular figures that should/shouldn’t be included in the main text?

  4. Do you feel a multiple linear regression using parameter values (rmax, Qc) and the covariates provided above is the best approach for a post-hoc analysis? Would creating a structural equation model (see Supplementary Figure 1) be a better approach?

  5. Is there anything else you feel I should consider examining moving forward?

Supplementary Figures
Supplementary Figure 1. Hypothesized structural equation model with primary covariates included in this project. Sources denoted numerically can be found listed at the end of this document.Supplementary Figure 1. Hypothesized structural equation model with primary covariates included in this project. Sources denoted numerically can be found listed at the end of this document.

Supplementary Figure 1. Hypothesized structural equation model with primary covariates included in this project. Sources denoted numerically can be found listed at the end of this document.

Supplementary Figure 2. Persistence curves for a site that performs well with the model structure described above (South Branch Potomac River) and a site that performs poorly (Reedy Creek). The critical discharge at which biomass is disturbed (c) is denoted by the dotted vertical line. The sensitivity of the persistence transition from presence to removal of biomass (s) is equal to the slope of the green persistence curve as it transitions from complete persistence (P = 1) to removal (P = 0) of biomass.

Supplementary Figure 3. Plots of median maximum growth rates (r~max~) versus (A) daily mean light availability, (B) mean daily summer water temperature, (C) watershed area, (D) stream order, (E) stream width, (F) latitude, (G) relative influence of dams, (H) mean nitrate concentrations, and (I) mean orthophosphate concentrations. Note, in panel G, a value of 0 indicates the certain influence of a dam while a value of 95 indicates a low or unlikely influence of a dam. Also, panels C, E, H, and I present the x-axis on a log scale.

Supplementary Figure 3. Plots of median maximum growth rates (rmax) versus (A) daily mean light availability, (B) mean daily summer water temperature, (C) watershed area, (D) stream order, (E) stream width, (F) latitude, (G) relative influence of dams, (H) mean nitrate concentrations, and (I) mean orthophosphate concentrations. Note, in panel G, a value of 0 indicates the certain influence of a dam while a value of 95 indicates a low or unlikely influence of a dam. Also, panels C, E, H, and I present the x-axis on a log scale.

Supplementary Figure 4. Estimates of the critical discharge at which biomass is disturbed (c) based on model fits at six sites using discharge normalized to either the ten-year flood at a given site based on the past 50 years of USGS records (Q10) or the maximum discharge in a given record (Qmax).

Supplementary Figure 5. Plots of the ratio of the median critical discharge (Q~c~) to the two-year flood (Q~2yr~) versus (A) stream order, (B) stream width, and (C) relative influence of dams. Note, in panel C, a value of 0 indicates the certain influence of a dam while a value of 95 indicates a low or unlikely influence of a dam. Also, panel B presents the x-axis on a log scale.

Supplementary Figure 5. Plots of the ratio of the median critical discharge (Qc) to the two-year flood (Q2yr) versus (A) stream order, (B) stream width, and (C) relative influence of dams. Note, in panel C, a value of 0 indicates the certain influence of a dam while a value of 95 indicates a low or unlikely influence of a dam. Also, panel B presents the x-axis on a log scale.